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Abstract 

The Aqua spacecraft will carry four single-axis gyros configured with three orthogonal axes and one skew 
axis. This redundancy presents a challenge for batch methods of on-orbit gyro calibration that use a 
spacecraft rotation model deterministically related to gyro data, in that sensor data can respond to at most 
three angular velocity components. When the number of gyros, N, is greater than 3, the 3xA/ matrix, G, 
that reduces the N gyro measurements to three body-frame angular-velocity components cannot be fully 
determined by such methods; there are many such matrices that produce essentially the same angular- 
velocity history. In such a case, spacecraft operators require information about the Nx3 gyro linear 
response matrix, R, that relates gyro outputs to the body-frame angular velocities causing them. This 
matrix provides sufficient information to determine multiple reduced-dimension G-matrices for use in case 
of failure or degradation of one or more gyros, as well as to determine an optimal 3xN G for the fully- 
functional configuration. 

A viable proposal is to apply a 3xA/ pre-filter matrix, F, to the N gyro outputs before carrying out a 
conventional gyro calibration procedure. The angular-velocity history emerging from conventional calib- 
ration may then be used as input data, together with the same gyro data that generated it, to fit the 
alignment, scale-factor, and bias parameters of each gyro axis in turn. A difficulty of such a proposal is 
the arbitrariness in the choice of F. Due to gyro noise, different pre-filter matrices produce different 
calibrations. This paper presents a method of choosing F that is based on optimizing gyro consistency in 
the limit of infinite weight on gyro data, as compared to sensor data. The choice of F is independent of a 
priori alignment and is based on the gyro data alone. The method is applicable to any N of three or more, 
but reduces to conventional batch-estimation methodologies when N=3. Results of computational 
comparison among calibration simulations using various choices of F will be presented for the Aqua gyro 
configuration with N=4. 


INTRODUCTION 

When the Aqua mission is launched later this year or early the next, it will carry with it one more gyroscope sensitive 
axis than is needed to fully determine the angular velocity. The purpose of the fourth gyro axis (for brevity, “fourth 
gyro”) is to provide redundancy in case of the failure of one of the others; current understanding is that it will be 
turned off awaiting need, in normal operation. Aqua can, however, be run with all four gyros turned on, connected 
to the control system, and reporting their measurements to the ground via telemetry. This presents the opportunity of 
using any of the gyros to check the accuracy of the other three. But it presents, as well, the challenge of resolving the 
discrepancies, due to finite measurement accuracy, that will inevitably arise. It will not, in general, be just one of the 
gyros that is responsible for the discrepancy, and there is no way uniquely to identify how the errors on the gyros 
have combined to produce an observed discrepancy. 

This would be all opportunity and no problem if there were accurate knowledge of the on-orbit angular velocities so 
that gyro calibration could be a process of fitting the calibration-curve slopes and intercepts to the observations of 
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gyro readouts and angular velocities, like calibrating a thermometer against a standard. Instead, the sensors against 
which one must compare the gyro readouts are attitude, not rate, sensors and relatively less accurate than the gyros 
(unless long time-intervals are used). The gyro calibration methods which have successful operational experience 
(Reference 1) at Goddard Space Flight Center (GSFC) all require that one be able to model the estimated angular 
velocity as a function of the gyro readouts and calibration parameters via an equation equivalent, in some form, to 

cb*=Gg-b (1) 

where ca is the estimate of the angular velocity vector in body-fixed coordinates, g is a column-vector of gyro 

readouts, G is a matrix of the appropriate size, and b is a bias vector (which may or may not have some slight pro- 
grammed time dependence). This expression is integrated to generate a family of models of some attitude observ- 
ables (such as attitude sensor data or slew offsets) depending on G, b and other model parameters such as initial 
attitudes. The calibration parameters can then be adjusted to fit the observables. At first thought, it appears neces- 
sary only to expand G from 3 columns to N. The problem is that, however many components g may have, the 
manifold of g , to the extent that the gyros respond only to rotation rates, is a space of, at most, three-dimensions 
mapped out by the variations of the true angular velocity. All G-matrices that produce equivalent results on that 
space are forever indistinguishable. 

Equation (1), in reducing the dimensions from N to three, implicitly resolves any internal disagreements in the gyros. 
The better the accuracy of the gyros, the less is the observable difference among competing models for doing that 
and the more ambiguous and subtle the choice among them becomes. In the end we find that the choice cannot be 
left purely to the retrospective comparison of attitude trajectories propagated from model angular velocities and 
compared to attitude sensor data. The extra dimensions of information in the gyro data, in other words, must not be 
discarded at the start by merely modeling the three-dimensional angular velocities. Yet we shall show that the 
existing three-axis calibration methods are still usable with the addition of a little superstructure to resolve the 
retrospective observation ambiguity by an imposed requirement of predictive optimization. 

MATHEMATICAL STATEMENT OF THE PROBLEM 

The N outputs at n times g, j (i=l...N, j=l...n) of a well behaved gyro assembly can be stacked into column N- 
vectors and put side by side in an Nxn matrix. The time sequence does not need to consist of continuous gyro 
readouts, but can include any number of separate time spans. The linear response of each gyro to the stimulus of 
angular velocity (as represented by n angular velocity 3-vectors lined up in the 3 xn matrix to) is described by the 
equation 

g = Roj + BU + v (2) 

where R, the gyro linear-response matrix, is Nx 3, B, the bias matrix, is Nxm, U is mxn, and V is Nxn gyro error, to 
be modeled as noise. [The form BU we have written for the biasing function requires some explanation. It is 
constructed to allow different constant biases for any of the separate gyro time spans. Each column of B is one of m 
independent bias ^/-vectors, each applicable to one or more of a set of m or more disjoint time intervals. The matrix 
U indicates which of the m biases corresponds to each of the n times. Matrix U is taken to be made up of n columns, 
each of m-1 zeros and a single 1. In each row of U, the l’s are grouped in clumps corresponding to the time inter- 
vals. More general forms of U might be used, but this corresponds to the type of time-dependent biasing function 
supported by the Batch Least-Squares IRU Calibration method (BICAL) (Reference (1)) and is convenient for this 
analysis.] The angular velocities, to, are written in an appropriately defined orthonormal body-fixed coordinate 
system, so that R and B may be approximately constant for strap-down gyros. 

One would then wish to solve this equation for an estimate, to*, of the body angular velocities in the form 

<o* =Gg-bU (3) 
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The angular velocity estimation error, to* - to, will be devoid of systematic variation with to, if and only if 

GR = I 3 (4) 

GB = b (5) 

When N= 3, these equations are sufficient to determine the gyro calibration parameters G and b for given gyro 
physics represented by R and B. When N>3, Equation (4) is insufficient to determine G for given R (or R for given 
G). R has many left inverses (and no right inverses). The existence of extra gyros is, however, a boon not a bane, 
and one can use them to reduce the effect of noise on to* by a factor something like V(3/7V). For equally noisy gyros, 
the optimal noise-reducing G (if R is known) is the well-known pseudoinverse of R 

G = pinv(R) = (R r R) 1 R T (6) 

which is one of the many matrices satisfying Equation (4). For sensors considered noisy relative to the gyros, 
however, one should not expect that any one of the left inverses of R would provide unambiguously best agreement 
with the sensor data. The residual disagreement will be dominated by the sensor error, whereas the differences 
provided by various left-inverses of R are of the order of gyro noise. 

With this background, we introduce a loss function for gyro calibration consisting of the sum of gyro and attitude- 
sensor terms: 

7 = / gy jR,B,<o)+7_ 0 >,a) (7) 

•f sensor measures, in a least-squares sense, the agreement of some complement of attitude-sensitive measurements or 
attitude estimates with the results of propagating the epoch attitudes a using the angular velocities to. It is the usual 
BICAL or Davenport-algorithm loss function, which we don’t need to specify further at this stage. The gyro loss 
function is 

7 gyro = ~Tr{ [Reo + BU - gf M[Rco + BU - g]W } (8) 

where W and M are symmetric positive semidefinite square matrices of dimension nxn and NxN giving the 
observation weights by time and by gyro-number, respectively. Although more general forms may be used, we shall 
take them to be diagonal matrices, corresponding to the assumption of lack of correlation among gyro measurements. 
Specifically, 

m=s t s 

S = diag([l/cr 1 ,.,.,1/crJ) (9) 

W = diag( [Af, , . . . , At n ] ) 

with the sigma’s being gyro noise specifications in radians/sec' /2 . This is a general least-square loss function for 
gyros, except for the assumptions that gyro observation noise has a factorable dependence on time and gyro 
instrument number, and that the time dependence is pure white rate-noise. 

This loss function depends on N(3+m)+3n+3n a independent variables (where n a is the number of attitudes a), which 
is very large when the number of gyro measurements, n, is large (20,000 or so in the Terra calibration effort). It 
might be doubted whether the minimum in such a large number of variables can be either well-defined or practically 
calculable. At any fixed R and B, however, the optimization with respect to each column of co is well-defined, even 
with only i gyro in play, and i sensor cannot make it less well-defined. Much of the information in J gyro will be consumed 
in that conditional minimization, but all of the information in J sensor remains available to constrain the remaining 
N(3+m)+3n a degrees of freedom of R, B, and a. We shall show proof, below, that the minimum of Equation (7) is 
well-defined in the limit of infinite weight on 7 gyro relative to 7 scnsor . So the minimum cannot be degenerate for all 
relative sensor/gyro weightings. As to the computability of the minimum, it can be stated, without proof here, that a 
change of variables can give the normal equations of the problem a sparsity structure permitting their solution by 
sequential methods (bearing a close resemblance to a sequential filter-smoother) with computational effort that grows 
only with the first power of n, for large n. Numerical stability and ultimate practicability of that solution are 
unproven at this time, and it is not the subject of this paper. 

In this paper, we shall accept the more modest goal of establishing contact with previous successful on-orbit gyro 
calibration methods by using the limit of infinite gyro weight to establish a form for a low-dimensionality family 
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{ co* } of approximate solutions for co. This limit of infinite gyro weight in the loss function must not be confused with 
the limit of zero gyro noise. It is distinct in that 

• it does require the loss-function to become non-optimal, of course, but this disadvantage is trumped by the 
fact that 

• it is realizable with real gyro data, which do not satisfy Equation (2) with zero noise. 

This limit corresponds to the assumptions made in the successful gyro calibration methods of Reference 1 for N=3. 
In that case, Equation (3) is the obvious solution to the infinite-weight limit, but 7 gyro drops out of the problem by 
achieving a minimum at exactly zero (taking R as the inverse of G) and so only 7 sensor requires explicit minimization. 
For N>3, that limit does not trivialize J gyto (unless the gyro data are noise-free simulated data). 

Our plan for the approximate minimization of the loss function of Equation (7), carried out in the succeeding 
sections, will then be as follows: 

• For N= 3, minimization of J gyro with respect to (0 implies: 

■ Equation (3) with G and b determined completely by Equations (4) and (5). 

■ G and b have the same information content as R and B, and may replace them as the independent variables 
of the problem, as in traditional treatments such as Reference 1 . 

■ The information content of 7 gyro is completely consumed in the minimization with respect to to, the 
minimum of J gym is exactly 0, and it disappears from the problem. This explains why 7 gyro need not appear 
explicitly when N- 3, as in Reference 1. 

• For N>3, unconditional minimization of 7 gyr0 with respect to to will be shown to require 

■ Equation (3) with G and b satisfying, but not completely determined by, Equations (4) and (5). 

■ A generalization of Equation (6) (see Equation (11) below) will be shown to be a sufficient (but only 
weakly necessary) condition for the minimization. 

■ R and B must be maintained as the independent variables and may not be replaced by G and b. 

• Some of the information content of J gyi0 survives the partial minimization and its retention is crucial to the 
solution of the problem. 

• The useful result that shall be established is the reduction of Equation (3), with A-column inputs g, to Equation 
(39), with 3-column inputs obtained linearly from g (Equation (40)) by multiplication with a 3 xA “filter” matrix 
F. The condition of minimizing the gyro loss function, subject to Equation (11), provides the value of F. 

• Any method applicable to the calibration of a 3-gyro complement (Reference 1) may then be used to adjust the 
parameters ^i 3 and b of Equation (39) to those that minimize the sensor loss-function, 7 sensor . That process 
provides the G-matrix as the product G 3 F. 

• If the gyros have sufficiently little noise, filter matrices F other than the one that minimizes can give almost 
equivalent angular-velocity sequences in Equation (39). G-matrices from such calibrations will not be optimal 
noise-reducing ones. Decent values of R and B will then result from a constrained minimization of J&y 0 with co 
replaced by these approximate angular velocities, since they differ little from the optimal ones. More nearly 
optimal G and b may be derived from these R and B determinations via Equation (11). 

• This process could be iterated and would be expected to converge to the results that would have been obtained 
immediately if the optimal F had been used from the beginning. 

SOLUTION IN THE LIMIT 

In the limit where unlimited faith is reposed in the gyros compared to the sensors, 7 sensor can only be minimized 
subject to the constraint that J gyro is already at its absolute minimum. (That minimum must be degenerate for 7 sensor to 
have any effect on the problem, but it is — see below). In particular, that assumption allows the angular velocities to 
be solved for in terms of R and B by using 0 = dJ gyio /d (0 =4> 

R r M(Roo* + BU - g) = 0 (10) 

The angular velocities solving this equation are given by Equation (3), provided that b is defined by Equation (5) and 
G by 

G =G pinv (R) = (R r MR) 1 R r M (11) 


178 



This slightly generalized (i.e. M-weighted) pseudoinverse form of G has the properties 

R = M _1 G 7 (GM~‘G r ) _1 (12) 

(RG) r MRG = MRG = (MRG) r (13) 

It should be emphasized that Equation (11) comes from a solution for to*, not from a solution for G. This solution is 
time-independent in that it can be established separately for each individual column of to, independently of the time- 
dependence and of the relative weighting (W). The G-matrix is only a surrogate for a family of co-solutions. In 
particular, for (nearly) noise-free gyro data when N>3, as previously discussed, many G-matrices (all those solving 
Equation (4)) will produce (nearly) the same to*. Nevertheless, this need not discourage one from using Equation 
(11) — it will always be able to produce the desired to’ even when that ability is not unique. Even when there might 
be other G-matrices that could produce marginally better fits to the retrospective data in the loss-function J, if one 
thinks one knows something about R, one is better off using the weighted pseudoinverse for G to reduce gyro data in 
the future . Perhaps Equation (11) should be viewed as a constraint on the problem, rather than something that can be 
determined from the data. 

Even apart from the limit of infinite-gyro-weight, Equation (3) may be introduced as an additional assumption to 
reduce the number of degrees of freedom in the loss function (by removing all the thousands of angular velocity 
components as independent variables). Substituting Equation (3) into (8) and (7): 

/(R,G,B)=7 g ,„(R,G,B)+y„„(G,b,<j) (14) 

V(R,G,B) = lTr{[(l N -RGXg-BUf M[(l N -RGXg-BU)]w } (15) 

Equation (1 1) is sufficient to guarantee zero for the gradients of 7 gyro with respect to G (and necessary, as well, unless 
the gyros are actually noiseless). 

Let us now minimize the loss function with respect to the m N-dimensional biases B: 

0 = ^gyroP B => 

(l N -RG) r M(l N -RGXg-BU)WU r =0 (16) 

(l N -RG)B = (l N -RGXg) (17) 

where 

(g)^gWU r (UWU r r 1 (18) 

is a matrix (of the same size as B) of the gyro measurements averaged over each biasing subinterval. Equation (17) 
supplies the information required in addition to Equation (5) to determine a B that minimizes the total loss-function 
from a b that minimizes only y sen sor- When that latter equation is substituted for GB in the left-hand side of the 
former equation we have 

B = (l N -RG)<g) + Rb (19) 

which satisfies Equation (5) thanks to Equation (4). Furthermore, substituting Equation (17) into Equation (15) 

V(R,G)=iTr{M(l N -RG]C ! (l n -RGf ) (20) 

where 

C g =<%W Sg T (21) 

5g = g- < g > U (22) 

Patently, ./ gyro is invariant under the transformation G— »AG together with R— >RA 1 (where A is an arbitrary 
nonsingular 3x3 matrix). This transformation preserves the mutual pseudoinverse relationship of R and G. 
Minimization of 7 gyro is degenerate over at least these nine degrees of freedom. Therefore the constraints provided 
by infinitely weighting 7 gyro relative to / sensor still allow us to minimize the latter using 
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G = G 3 F (23) 

where G 3 is a 3x3 matrix which may be arbitrarily varied along with b to optimize the agreement with sensor 
observations, while the 3 xN matrix F is responsible for maintaining the optimization of the gyro loss function. Note 
that Fg is a set of 3-dimensional gyro-derived data. The minimization of 7 sen sor with respect to G 3 and b is precisely 
the same problem faced in the calibration of spacecraft with 3 gyros and can be treated with any of the existing 
methods, for example, BICAL. The new problem is now to find the F that optimizes 7 gyro . 

As a preliminary to solving that problem, let us digress upon the meaning and properties of the matrix C g . It is 
proportional to the sample covariance matrix (over the sample of gyro observations) of g, with means subtracted on a 
subinterval-dependent basis. Since the g are constrained to obey an equation of the form of Equation (3), with 
angular velocities of three dimensions only, one should expect this matrix to have only 3 eigenvalues substantially 
differing from zero. Consider a set of gyro data from a truth model given by Equation (2) with R=R t and w=to,. 
Then 

SCS r /« = S (R t &> t +£v)W(R t <5co, +Sv) T S T /n 

8 (24) 

-S R t ((5o3 t W(5(o t /n)R t r S r + I N 

Assuming that the set of angular velocities in the calibration dataset is adequately three-dimensional, the first matrix 
in the expression above is of rank three having N - 3 null eigenvalues and 3 of order (8o)/foAf l/2 ]) 2 . We assume, 
as we must if we are to expect any gyro calibration to yield good results, that the spacecraft has been maneuvered 
about all axes at statistically significant rates, so that the three nonzero eigenvalues are large compared with 1 , the 

eigenvalue of the noise term approximating I v . The addition of the second (identity) term ensures that SC g S 7 is 
nonsingular with 3 large eigenvalues and N- 3 much smaller. 

Varying 7 gyro (Equation (20)) with respect to G at fixed R (not yet assuming Equation (11)) 

o=ay g>ro /3Gi R => 

R r M (l N - RG) C g = 0 (25) 

Equation (1 1) is implied by this only to the extent that C g is nonsingular. As we have seen, C g is nearly singular in 
N- 3 dimensions, and can in fact be singular if there is no gyro noise, such as with error-free simulated gyro data. 
This simply reminds us that the choice of a particular left-inverse of R for G is of secondary importance to the 
necessity for GR to be I 3 . The choice of G pinv (R) is motivated by noise reduction and is not important if the noise is 
actually absent. Nevertheless, even when C g is singular, this choice of G still does satisfy Equation (25) and 
minimize the loss function — it simply no longer does so uniquely. Of course, it only does so exactly in the limit of 
infinite gyro weight. Since some of the gradients of 7 gyro , i.e. those in the direction of the small-eigenvalue 
eigenvectors of C g , would be small even for other Gs, the influence of the gradients of 7 sensor cannot necessarily be 
neglected. We will have gained little if we consider that we have added 7 gyro (R,G,B) to 7 sensor (G,b,a) along 

with the “new” variables R and B; the power of 7 gyro is largely used up in defining R and B as functions of the old 
variables G and b — a useful thing, since R can imply G-matrices to use for different weightings, M, and for 
degraded 3-gyro configurations, but it won’t make the minimum much more unique. We have to remind ourselves, 
however, that we now view G and b as but proxies for co', and that it is quite reasonable to impose constraints on G 
so that it does not import nearly irrelevant new destabilizing degrees of freedom into the problem. 

Now let us finally work out the constraints on R corresponding to minimizing 7 gyro with respect to R at the 
corresponding minima of G and B. Substituting Equation (13) into (20). 

7»»( R ) = y Tr { M [!n ~ RG l C g } M G = G pl „(R) (26) 

Minimizing this is equivalent to maximizing the benefit function 

A„=Tr{MRGC ! } (27) 

In order to find the constrained (by Equation (11)) maximum of 0 , we now consider the “economy-size” singular 
value decomposition (SVD) of SR: 
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SR = V R D R U R (28) 

where 

U R = 3x3 unitary 

D r =3x 3 diagonal, nonsingular (29) 

V R = N x 3 , 3 columns of an N x N unitary matrix, V R T V R = I 3 V R V R r ^ I iV 
Using this SVD we find, after some algebra, 

G = U r D r _ 1 V r S (30) 


MRG = S' V R D R U' R U R D R _1 V 7 S 


r V R D R U ri 

r v R v 7 

( 3 


=s'y r \ t r s 


= S T 




2> R y R ,. 

V i=1 

where v Ri are the orthonormal columns of V R . Substituting this into (27) 

5 gyro=X V R! SC g Sry R<- 


(31) 


(32) 


The maximization of this, subject to the orthonormality constraints on the columns of V R , has the general solution 


v R = v c u 3 ' 


where U 3 is an arbitrary 3x3 unitary matrix and 

V c =[v 


Cl T C2 


v cs] 


(33) 


(34) 


is composed of columns which are the orthonormal eigenvectors corresponding to the largest three eigenvalues of 

SC g S r : 


SC g S r V c = V C D C (35) 

D c is a 3x3 diagonal matrix with entries » n. It is not necessary that eigenvalues be nondegenerate; these eigen- 
vectors do not need to be cleanly distinguished from each other. They do need to be clearly distinguished from the 
eigenvectors of the N - 3 smaller eigenvalues. What is important is the subspace spanned by the three eigenvectors; 
any mixture of the top three eigenvectors may be absorbed into the unitary matrix U 3 . 


Then the Gs that maximize f? gyro are of the form 


where 


G = U r D r - ! U 3 V c 7 S 

=g 3 f 

(36) 

g 3 = u r d r - , u 3 a- 1 

(37) 

F = AV/S 

(38) 


for an arbitrary nonsingular 3x3 matrix A (which will be used, below, to match a priori conditions). The expression 
for G 3 A, above, is the SVD of a perfectly general 3x3 matrix, there being no constraints on U R , D R , U 3 other than 

those defining the SVD. G 3 is not constrained or restricted by the minimization of the gyro loss function. It is the 
3xN matrix F that bears the full responsibility for that minimization. 


We are therefore free to use G 3 and b to minimize the sensor loss function with the equation 

to* = G 3 to R -bU (39) 

where 
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Or = Fg (40) 

is a three-dimensional reduction of the gyro measurements. It is just this problem that the existing 3-dimensional 
BICAL algorithm solves, as well as other conventional 3-dimensional gyro calibration methods. The resulting of 
could be substituted into Equation (2) to provide data for a linear fit of that equation to determine the m degrees of 
freedom of B and 3(N-3) of R and that are undetermined by Equations (4) and (5). But since this linear fit is already 
implicit in the gyro loss function, we get equivalent results by minimizing Equation (20) with respect to R. Since we 
satisfy Equation (25) with Equation (1 1), we do not need to worry that G varies implicitly with R in now finding the 
global minimum of J gyr0 by imposing 0 = dJ gyi0 /d R = a/ gyro/ aRl c => 

S(l N -RG)C g G r =0 (41) 

with the solution 

R = C g G r (GC g G r ) 1 = RG 3 1 (42) 

R = C g F r (FC g F r ) _1 (43) 

The usefulness of the above equation is not restricted to the infinite-gyro- weight limit, nor to F being optimal, i.e., as 
given by Equation (38). One might, for instance, define a set of conservative calibration strategies by assuming 
particular fixed values for F. For example, if F were taken to be the first 3 rows of I 4 , the sensor loss-function fit 
would be a conventional calibration of the first 3 gyros. This would provide the first 3 rows of the R-matrix as G 3 '', 
and the first three rows of B as Gf'b, but the fourth rows would have to be found by comparing the fourth gyro’s 
measurements to the angular velocities of. Equation (42), with the help of Equation (19), performs that task, and it 
is easy to show, with those equations, that FR (in this example, the first 3 rows of R) is G 3 ‘ 1 and FB is G 3 ''b. 

We now have another way to derive the optimality condition for F (Equation (38)) as being the condition of 
consistency between Equation (11) (the requirement of best noise reduction by G for known R) and Equation (42) 
(the requirement of best fit of the gyro data to (o’ for known G). For, if both these conditions are true, then the 
partials of J %yi0 with respect to both G (at fixed R) and R (at fixed G) are zero, as shown above, J gyro must be at a 
minimum and Equation (38) follows. And indeed, Equation (42) can be shown consistent with Equation (11), and 
the proof uses Equation (35) essentially. 

If F is optimal and thus of the form given by Equation (38) we have, substituting that equation into Equation (43), 

SR = SC g S T V c A T (AV c r SC g S r V c A T r 1 

= v c d c (v/v c d c )- 1 a - 1 

= V C A-' 

And then 

Gptov(R) - ((sr( sR^sity s 

= |v c A-') r V c A- , ) r, (v c A- , ) r S 

= AV c r S = F 

G pinv (R) = G pmv (RG 3 -') = G 3 G pinv (R) = G 3 F = G 

The set of conservative calibration procedures that pre-defines F, not taking advantage of optimal noise reduction in 
the 3-dimensional 7 se nsor calibration problem, might proceed as follows. Minimize 7 se „ sor as a function of G 3 and b. 
Then fit R and B to the estimated of implied by Equation (39) by using Equations (42) and (19). With this R in 
hand, one would want to support the mission’s future gyro data using optimal noise reduction, so one would define 

G new = G 3 G pinv (R) = G 3 F new (46) 

for this purpose. One might think that there could be some gain in iterating this procedure further, each time 
replacing F according to the above and re-minimizing 7 sen SO r with respect to G 3 and b. This would be burdensome, 


( 44 ) 


( 45 ) 
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since the minimization of 7 se nsor is by far the most computationally intensive part of the process. Such a procedure 
would, however, probably converge, and, if it did, it would have to be to an optimal F according to Equation (38). 
Using the optimal F, determined in advance from eigen-analysis of the gyro-only data matrix C g , permits an 
equivalent answer to be achieved with a single J sc „ sor minimization. 

An a priori value for G 3 is needed to initialize the 3-dimensional calibration (by a method such as BICAL or the 
Davenport algorithm) using the sensor loss-function. This may be obtained conveniently using Equation (4) if one 
has an indication of a nominal R-matrix, R°. An a priori G 3 of G 3 °=I 3 can be used if one chooses 

A = (V/SR V (47) 

(or use this expression for G 3 ° with A=I 3 ). 

The /V-row matrices R and B provide full information on all the gyros’ response to spacecraft rotation. Knowledge 
of them permits us to make recommendations for versions of the G-matrix (and three-dimensional biases b) to be 
used when any one of the gyros must be down-weighted because of increased noise or even totally excluded from 
use. One simply uses Equations (11) and (5) with the same (calibrated) R and B but with a new weight matrix M new . 
If a row and column of M new is set to be zero, the corresponding gyro will have a null column in the resulting 
weighted pseudoinverse, so that that gyro’s data are removed from further use in calculating angular velocities. A re- 
calibration is not necessary (unless, of course, a partially degraded gyro is thought to have changed its biases or 
sensitivities). 

RELATION TO THE SOLUTION FOR OPTIMAL WEIGHTING 

The limit of infinite weight on 7 gyro versus 7 se nS0 r has been used in the above section. Strong arguments can be made 
that the minimum of Equation (14) is little dependent upon the relative weighting of the two terms. (Little, however, 
can be said about the use of the infinite weight limit to get from Equation (7) to that point.) J sem0 r has only small 
sensitivity, vanishing with gyro noise, to variations dG of G for which (dG)R t =0. 7 gyro , on the other hand, has full 
sensitivity to variation in those directions proportional to the inverse square of the gyro noise parameter. 
Furthermore, the 12-dimensionally degenerate minimum of ./ gyr0 can easily be attained using only such variations at 
little cost to ./ S ensor- Under such circumstances, the relative weighting of the two terms has marginal effect on the 
solution. A fully mathematical demonstration of this is possible, but requires a large expansion of the notation and 
will not be presented. 

SIMULATED NUMERICAL EVALUATION 

We have tried these methods on simulated data corresponding to the Aqua spacecraft, with all four gyros operating 
and reporting data, while undergoing a calibration maneuver sequence like that executed by the Terra mission on 
April 13, 2000 (Reference 1). The eight calibration slews started with a positive (x-axis) 10-degree roll followed by 
a roll back to nominal. This sequence was repeated in the negative roll direction. There were similar 30-degree 
positive and negative yaw (z-axis) offset maneuvers. Terra, like Aqua, is an Earth-oriented spacecraft and at all 
times it maintained its 1 rotation per orbit attitude motion about orbit normal. Because of this constant rotation, the 
offset attitudes resulted in changes of the pitch rates by the cosine of the offset angle and therefore provided rather 
weak, but realistic, pitch (y-axis) scale-factor observability. This sequence was approximated by simulated rate data 
with the nominal body y-rate of -0.06 deg/sec. Slews were performed at 0. 1 deg/sec, offset attitudes were held for 
1000 sec with intermediate nominal holds for 5000 sec, and the whole sequence was begun and ended with a 1000- 
sec hold at nominal rate. The entire sequence required 22600 sec. The angular velocities corresponding to this 
scenario were recorded every second (and each second had constant rates). 

This angular-velocity history was integrated exactly to provide a truth model of the attitude history. Attitude 
quaternion “measurements” every 8 seconds stood in for sensor data, with pseudorandom Gaussian noise standard 
deviation of 10 arcsec per axis. Gyro data were formed from Equation (2) with a truth R-matrix, R t corresponding to 
the nominal for the Aqua mission: 
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This is an orthonormal triple of gyros supplemented by a fourth gyro sensitive to the relative [-1-1 -Indirection, 
with the whole assembly rotated to put the fourth axis in the z (yaw) direction. The angle between gyro axis 4 and 
the others is 125.26 deg, rather than the 109.47 deg that it would be for a regular tetrahedral configuration. 
Pseudorandom Gaussian noise impressed on the gyro data included white rate noise at three different levels (4, 57, 
and 229 xl() 6 deg/sec 1/2 ) and bias random walk (6xl0‘ 9 deg/sec 3/2 ). The lower level of gyro noise realistically 
corresponds to the specification for a NASA standard gyro of the UARS/EUVE era, the Dry-Rotor Inertial Reference 
Unit, Model II. The other two levels are moderately to highly pessimistic, and were included in order to exacerbate 
the differences among different calibration methods for illustrative purposes. All three simulations used the same 
random-number samples, scaled differently. 


For each gyro noise level, four different calibration methods were exercised. All used BICAL to minimize 7 sensor 
with respect to G 3 and b. 

• The “optimal” method allowed F to vary and find its optimum (Equation (38)) with G=G pinv (R) in the 
infinite-gyro- weight limit. 

• The “suboptimal” method fixed F at the pseudoinverse of a nominal R° close to R t . (This R° takes the 
place of pre-launch calibration assumptions and has to differ from R t . Gaussian pseudorandom numbers of 
standard deviation 0.002 were added to R t to form R° with scale-factor errors ranging from -1.6 to +5.1 
parts per thousand and alignment errors ranging from 2.1 to 3.9 milliradians (0.12 to 0.23 degrees).) 
Calibrated R-matrices were then calculated, in all cases, with Equation (42) (linear fit to the estimated 
angular velocities using G=G 3 F). 

• The third calibration method (called “no-4”) set F at the left-inverse of R° that ignores the fourth gyro axis, 
doing a three-axis calibration of the remaining orthonormal triple of gyro axes, before fitting the fourth gyro 
to its data. 

• The last method (called “no-3”) used a fixed F analogous to the “no-4” method, but with gyro axis 3 rather 
than 4 initially discarded. 

The differences among these methods would be nil except for gyro error. 

Table 1 shows the results for the maximum R-matrix element error in each simulation. 


Table 1. Summary of Simulated Calibration R-Matrix Errors 


R-Matrix Errors 

Rate Noise 
(microdeg/s 1/2 ) 

Calibration Method 

Optimal 

Suboptimal 


no - 3 

4 

1.4e-4 

1.4e-4 

1.6e-4 

1.2e-4 

57 

4.6e-4 

4.6e-4 

5.1 e-4 

2.0e-3 

229 

1.5e-3 

1.5e-3 

1 ,6e-3 

2.6e-2 


The optimal and suboptimal methods offer little difference; their R-matrices never differ by more than 12 parts per 
million. This is not unexpected, as the difference between these two is proportional to the product of gyro noise and 
a priori alignment error. The method that initially discards gyro axis number 4 also differs little from these two; that 
is perhaps more surprising because this foregoes the possibility of a 50% gyro noise variance reduction in the y-axis 
during the BICAL fit. The calibration method discarding gyro axis 3 performs poorly compared to the others, except 
at the realistic , low gyro noise level. All methods based on this paper perform quite well (0.01%) at realistic levels 
of noise. 
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Gyro rate noise dominates the errors at the two higher levels of noise, and, in view of that the poor performance of 
the “no-3” method is not so surprising. Consider the angular- velocity noise error 1 -sigma ellipsoids in three 
dimensions, assuming unit uncorrelated rate noise in each of the four gyros. For the nominal Aqua configuration, 
this is a spheroid, oblate in the z-axis, with semiminor axis of 0.71 (l/v/2) and unit semimajor axes. If axis four is 
excluded, the figure becomes a unit sphere. If one of the three orthogonal axes is excluded, on the other hand, the 
resulting canted triaxial ellipsoid has semiminor axis 0.74 and semimajor axis 2.3. This difference depending on 
which gyro is discarded is a bit surprising if one has been thinking that the difference between the Aqua 
configuration and a regular tetrahedral configuration is minor. For a regular tetrahedral configuration, the nominal 
error surface is a sphere of radius 0.87. Loss of any gyro only stretches that surface by a factor of 2 along the lost 
axis, only, to semimajor axis 1.7 (V 3). 

CONCLUSIONS 

Supernumerary gyro axes present a challenge for batch methods of on-orbit gyro calibration that use a spacecraft 
rotation model deterministically related to gyro data, but one that they can meet without fundamental changes. The 
3 xN matrix, G, that reduces N measurements to three body-frame angular-velocity components cannot be directly 
determined by such methods because the needed information is filtered out from the beginning by such an approach. 
But, by re-inserting the full suite of gyro measurements as observations in a least-squares loss function, the Nx 3 gyro 
linear response matrix, R, and a full N biases, can be determined from the results of 3x3 calibration. 

A variable element in this approach is the 3x/Vpre-filter matrix, F, applied to the N gyro outputs to give a 3-axis data 
stream to which ordinary gyro calibration is applied. The results of our simulations with N = 4 show that, for 
reasonable levels of gyro error, any well-founded choice performs adequately. The method of letting the gyro data 
determine an optimal F, developed in this work, certainly does well, but does not perform noticeably better than the 
use of the a priori nominally optimal F. Even without taking any advantage at all of the noise-reduction possible 
using the extra axes, such as by performing the 3x3 calibration with only 3 of the gyros, good results are usually 
achieved if the gyro to be discarded is chosen intelligently. This is not altogether surprising since, with only one 
extra axis, the amount of noise reduction is at most a factor of V 0.75. Use of the optimal choice of F would probably 
be of greater consequence for larger amounts of redundancy, although even for N=6 the potential for gyro noise 
reduction is only 30%. 

Certain facets of the method remain to be determined. Covariance analysis of the method is incomplete and will not 
be presented at this time. Since the infinite-gyro-weight limit for which exact results have been derived is a non- 
optimal smoother, there are some formal complications. A related issue is that of how to properly insert a priori 
information into the loss functions. Such information complicates the solution since it induces correlations between 
the 3x3 (G 3 ) and (N- 3)x3 (F) degrees of freedom that are absent from the on-orbit information, at least when the 
time-dependence (W) and gyro-axis dependence (M) of the gyro noise amplitude are factorable as we have assumed 
herein. 
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